function    A303()
format long;
% Markov链Monte Carlo（Markov chain Monte Carlo，MCMC）通过构造Markov链极限平稳分布来模拟计算积分。
% 遍历的Markov链的极限分布是唯一的平稳分布

cauchy = @(theta) 1 ./(1+theta.^2);                                          % 定义函数cauchy(theta)

% Metropolis算法
% 初始化Metropolis采样器
T = 10000;                                                                   % 迭代的最大值幅值
sigma = 1;                                                                   % 为正态密度分布标准偏差幅值
thetamin = -30; thetamax = 30;                                               % 定义初始值
theta = zeros(1, T);                                                         % 样本空间的初始化
seed = 1; rand('state', seed); randn('state', seed);                         % 随机种子赋值 
theta(1) = unifrnd(thetamin, thetamax);                                      % 产生初始值
% 开始采样
t=1;
while t < T                                                                  % 有T样本时，迭代
    t=t+1;
    % 使用可能的正态分布为theta设置一个新值
    theta_star = normrnd(theta(t-1), sigma);
    % 计算验收比率
    alpha = min([1 cauchy(theta_star)/cauchy(theta(t-1))]);
    % 归一化偏差[0 1]
    u = rand;
    % 是否接收该建议？
    if u < alpha
        theta(t)= theta_star;                                                % 若满足条件，样本赋值新状态
    else
        theta(t)= theta(t-1);                                                % 若不满足条件，采用前一个状态
    end
end
% Display histogram of our samples 
subplot(1, 3, 1);
nbins = 200;
thetabins = linspace(thetamin, thetamax, nbins);
counts = hist(theta, thetabins);                                             % 样本直方图
bar(thetabins, counts/sum(counts), 'k');                                     % 归一化直方图
hold on
xlim([thetamin thetamax]);
xlabel('\theta');
ylabel('p(\theta)');
title('Metropolis算法');
legend('样本直方图')
hold off

% Metropolis-Hasting算法
% 目标（真实）pdf为p(x)，这是不能采样的数据
% 采样（建议）pdf为q(x*|x)=N(x, 10)
X(1)= 0;
N = 1e4;
p = @(x) 0.3*exp(-0.2*x.^2)+0.7*exp(-0.2*(x-10).^2);
% 目标p(x)的可视化
dx = 0.5;
xx = -10:dx:20;
fp = p(xx);
% MH算法
sig = (10);
% sig = (1);
for i = 1:N-1
    u = rand;
    x = X(i);
    xs = normrnd(x, sig);                                                    % 新示例xs基于采样PDF中存在的X。
    pxs = p(xs);                                                             % 计算目标分布在建议样本处的概率密度
    px = p(x);                                                               % 计算目标分布在当前样本处的概率密度
    qxs = normpdf(xs, x, sig); qx = normpdf(x, xs, sig);                     % 产生p, q
    % if u < min(1, pxs*qx/(px*qxs))                                         % case 1: 伪码
    if u < min(1, pxs/(px))                                                  % case 2: Metropolis算法
    % if u < min(1, pxs/qxs/(px/qx))                                         % case 3: 独立采样器
        X(i+1) = xs;
    else
        X(i+1) = x;
    end
end
% 比较模拟结果的PDF与真实PDF
N0 = 1; 
nb = histc(X(N0+1:N), xx);
subplot(1, 3, 2)
bar(xx+dx/2, nb/(N-N0)/dx);                                                  % 样本可视化
hold on
A = sum(fp)*dx;                                                              % 归一化目标分布
plot(xx, fp/A, 'r')                                                          % 模拟结果的可视化
% 比较cdf和真实cdf。
F1(1) = 0; 
F2(1) = 0;
for i = 2:length(xx)
    F1(i) = F1(i-1)+nb(i)/(N-N0);                                            % 样本分布的 CDF
    F2(i) = F2(i-1)+fp(i)*dx/A;                                              % 目标分布的 CDF
end
plot(xx, [F1' F2'])
title('Metropolis-Hasting算法')
legend('样本直方图', '模拟pdf', '真实cdf', '模拟cdf')
max(F1-F2)                                                                   % 准确性衡量标准
hold off

% Gibbs算法
rng('default')
num_samples = 5000;
num_dims = 2;
mu = [0, 0];
rho(1) = 0.8; rho(2) = 0.8;
prop_sigma = 1;
minn = [-3, -3]; maxx=[3, 3];
x = zeros(num_samples, num_dims);
x(1, 1) = unifrnd(minn(1), maxx(1));
x(1, 2) = unifrnd(minn(2), maxx(2));
t = 1;
dims = 1:num_dims;
while t < num_samples
t = t+1;
T = [t-1, t];                                                                % 时刻信息的维护，上一时刻T(1)，下一时刻T(2)
for iD = 1:num_dims
    not_idx = (dims ~= iD);
    mu_cond = mu(iD)+rho(iD)*(x(T(iD), not_idx)-mu(not_idx));
    sigma_cond = sqrt(1-rho(iD)^2);
    x(t, iD) = normrnd(mu_cond, sigma_cond);
end
end
subplot(1, 3, 3)
h1 = scatter(x(:,1), x(:,2),'r.');
hold on
for t = 1:50
    plot([x(t,1), x(t+1,1)], [x(t,2), x(t,2)], 'k-');                        % x轴方向移动
    plot([x(t+1,1), x(t+1,1)], [x(t,2), x(t+1,2)], 'k-');                    % y轴方向移动
    h2 = plot(x(t+1, 1), x(t+1,2), 'ko');
end
h3 = scatter(x(1,1), x(1,2), 'go', 'Linewidth', 3);
legend([h1, h2, h3], {'Samples', '1st 50 samples', 'x(t=0)'}, 'location', 'northwest');
hold off;
xlabel('x_1');
ylabel('x_2');
axis square